MEDB 5501, Module03

---
title: "Normal probabilities and quantiles"
format: 
  html:
    embed-resources: true
editor: source
execute: 
  error: true
---

This program displays the standard normal curve along with various probabilities and quantiles. It was written by Steve Simon on 2024-09-01 and is placed in the public domain.

## Load the tidyverse library

```{r setup}
#| message: false
#| warning: false
library(tidyverse)
```

## Using R to draw the standard normal curve

use seq to calculate 100 evenly spaced values between -4 and +4 and dnorm to compute the bell curve at each point. Use geom_polygon to paint the area surrounded by the bell curve.

```{r standard-normal}
x <- seq(-4, 4, length=100)
y <- dnorm(x)
data.frame(x, y) |>
  ggplot(aes(x, y)) + 
    geom_polygon(fill="white", color="black") -> normal_curve
normal_curve
```

## P[Z < 1.5]

Use geom_vline and geom_label to draw a vertical reference line and add text to the normal curve. The pnorm function computes the standard normal probability.

```{r prob-1}
a <- 1.5
normal_curve +
  geom_vline(xintercept=a) +
  geom_label(x=a, y=0.4, label=a) +
  geom_label(x=a-0.5, y=0, label="Area = ?")
pnorm(1.5)
```

The probability that a standard normal variable is less than 1.5 is about 93%.

## P[Z < -0.5]

```{r prob-2}
a <- -0.5
normal_curve +
  geom_vline(xintercept=a) +
  geom_label(x=a, y=0.4, label=a) +
  geom_label(x=a-0.5, y=0, label="Area = ?")
pnorm(-0.5)
```

The probability that a standard normal variable is less than -0.5 is about 31%.

## P[Z > 1]

When you are calculating the probability on the right (probability greater than some number), use 1-pnorm.

```{r prob-3}
a <- 1
normal_curve +
  geom_vline(xintercept=a) +
  geom_label(x=a, y=0.4, label=a) +
  geom_label(x=a+0.5, y=0, label="Area = ?")
1- pnorm(1)
```

The probability that a standard normal variable is larger than 1 is about 16%.

## P[Z > -2]

```{r prob-4}
a <- -2
normal_curve +
  geom_vline(xintercept=a) +
  geom_label(x=a, y=0.4, label=a) +
  geom_label(x=a+0.5, y=0, label="Area = ?")
1- pnorm(-2)
```

The probability that a standard normal variable is larger than -2 is about 98%.

## P[-2.5 < Z < 2.5]

When you are calculating the probability between two values, compute pnorm of the larger value minus pnorm of the smaller value.

```{r prob-5}
a <- 2.5
normal_curve +
  geom_vline(xintercept=-a) +
  geom_vline(xintercept= a) +
  geom_label(x=-a, y=0.4, label=-a) +
  geom_label(x= a, y=0.4, label= a) +
  geom_label(x=0, y=0, label="Area = ?")
pnorm(2.5)- pnorm(-2.5)
```

The probability that a standard normal variable is between -2.5 and 2.5 is aobut 99%.

## P[-0.5 < Z < 0.5]

```{r prob-6}
a <- 0.5
normal_curve +
  geom_vline(xintercept=-a) +
  geom_vline(xintercept= a) +
  geom_label(x=-a, y=0.4, label=-a) +
  geom_label(x= a, y=0.4, label= a) +
  geom_label(x=0, y=0, label="Area = ?")
pnorm(0.5) - pnorm(-0.5)
```

The probability that a standard normal variable is between -0.5 and 0.5 is about 38%.

## 25th percentile of a standard normal

Use qnorm to calculate quantiles of the standard normal distribution.

```{r quantile-1}
p <- 0.25
a <- qnorm(p)
normal_curve +
  geom_vline(xintercept=a) +
  geom_label(x=a, y=0.4, label="Quantile = ?") +
  geom_label(x=a-0.5, y=0, label=p)
qnorm(0.25)
```

The 25th percentile of a standard normal distribution is -0.67.

## 90th percentile of a standard normal

```{r quantile-2}
p <- 0.9
a <- qnorm(p)
normal_curve +
  geom_vline(xintercept=a) +
  geom_label(x=a, y=0.4, label="Quantile = ?") +
  geom_label(x=a-0.5, y=0, label=p)
qnorm(0.9)
```

The 90th percentile of a standard normal distribution is 1.28.

   x r
1  7 4
2  3 2
3 23 9
4  2 1
5  5 3
6 13 6
7 11 5
8 17 7
9 19 8

   x r pctile
1  7 4    0.4
2  3 2    0.2
3 23 9    0.9
4  2 1    0.1
5  5 3    0.3
6 13 6    0.6
7 11 5    0.5
8 17 7    0.7
9 19 8    0.8

   x r pctile          z
1  7 4    0.4 -0.2533471
2  3 2    0.2 -0.8416212
3 23 9    0.9  1.2815516
4  2 1    0.1 -1.2815516
5  5 3    0.3 -0.5244005
6 13 6    0.6  0.2533471
7 11 5    0.5  0.0000000
8 17 7    0.7  0.5244005
9 19 8    0.8  0.8416212

---
data_dictionary: fev (.csv, sas7bdat, .sav, .txt)

copyright: >
  The author of the jse article holds the copyright, but does not list
  conditions under which it can be used. Individual use for educational
  purposes is probably permitted under the Fair Use provisions of
  U.S. Copyright laws.

description: >
  Forced Expiratory Volume (FEV) in children. The data was collected 
  in Boston in the 1970s.

additional_description:
  https://jse.amstat.org/v13n2/datasets.kahn.html

download_url:
  https://www.amstat.org/publications/jse/datasets/fev.dat.txt

format:
  csv: comma delimited
  sas7bdat: proprietary (SAS)
  sav: proprietary (SPSS)
  txt: fixed width
  

varnames: 
  not included
  

missing_value_code:
  not needed
  

size:   
  rows: 654
  columns: 5

vars:
  age:
    scale: ratio
    range: positive integer
    unit: years
    

  fev:
    label: Forced Expiratory Volume
    scale: ratio
    range: positive real
    unit: liters
    

  ht:
    label: Height
    scale: positive real
    unit: inches

  sex:
    value:
      F: Female
      M: Male

  smoke:
    value:
      N: Nonsmoker
      Y: Smoker
---

---
title: "Analysis of fev data"
format: 
  html:
    embed-resources: true
editor: source
execute: 
  error: true
---

This program assesses the normality of variables in a study of pulmonary function in children. There is a [data dictionary][dd] that provides more details about the data. The program was written by Steve Simon on 2024-09-02 and is placed in the public domain.

[dd]: https://github.com/pmean/datasets/blob/master/fev.yaml

## Libraries

The tidyverse library is the only one you need for  this program.

```{r setup}
#| message: false
#| warning: false
library(tidyverse)
```

## List variable names

Since the variable names are not listed in the data file itself, you need to list them here.

```{r names}
fev_names <- c(
    "age",
    "fev",
    "ht",
    "sex",
    "smoke")
```

## Reading the data

Here is the code to read the data and show a glimpse. 

```{r read}
fev <- read_csv(
  file="../data/fev.csv",
  col_names=fev_names,
  col_types="nnncc")
glimpse(fev)
```

## Calculate mean and standard deviation for fev

To orient yourself to the data, calculate a few descriptive statistics.

```{r descriptive-fev}
fev |>
  summarize(
    fev_mean=mean(fev),
    fev_stdv=sd(fev))
```

The mean fev is 2.6 liters and the standard deviation is 0.84 liters. I am not an expert on pulmonary function, but these values do not appear to be unusually large or small.

## Histogram for fev, wide bars

The geom_histogram function draws a histogram. You should specify values for color (which is the outline of individual bars) and fill (which is what is inside the bars). Also be careful with your choice of binwidth. Don't rely on the default choice.

```{r histogram-fev-wide}
ggplot(data=fev, aes(x=fev)) +
  geom_histogram(
    binwidth=0.5,
    color="black",
    fill="white")
```

See below for interpretation

## Histogram for fev, narrow bars

```{r histogram-fev-narrow}
ggplot(data=fev, aes(x=fev)) +
  geom_histogram(
    binwidth=0.1,
    color="black",
    fill="white")
```

Although some may interpret these histograms as showing a slight skewness, I would interpret them as being approximately normal.

## Normal probability plot for fev

The qqnorm function produces a normal probability plot. The default option for most plots is landscape orientation (the width is larger than the height). The q-q plot, however, looks best if figure width and height are equal.

```{r qqplot-fev}
#| fig-width: 5
#| fig.height: 5
qqnorm(fev$fev)
```

The normal probability plot is reasonably close to a straight line, indicating that the data comes reasonably close to following a normal distribution.